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ON THE FORMATION OF HOT JUPITERS IN STELLAR BINARIES 
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ABSTRACT 

We study the production of Hot Jupiters (HJs) in stellar binaries. We show that the "eccentric 
Kozai-Lidov" (EKL) mechanism can play a key role in the dynamical evolution of a star-planet-star 
triple system. We run a large set of Monte Carlo simulations including the secular evolution of the 
orbits, general relativistic precession, and tides, and we determine the semi-major axis, eccentricity, 
inclination and spin-orbit angle distributions of the HJs that are produced. We explore the effect 
of different tidal friction parameters on the results. We find that the efficiency of forming HJs when 
taking the EKL mechanism into account is higher then previously estimated. Accounting for the 
frequency of stellar binaries, we find that this production mechanism can account for about 30% 
of the observed HJ population. Current observations of spin-orbit angles are consistent with this 
mechanism producing ~ 30% of all HJs, and up to 100% of the misaligned systems. Based on the 
properties of binaries without a HJ in our simulations, we predict the existence of many Jupiter-like 
planets with moderately eccentric and inclined orbits and semi-major axes of several AU. 



1. INTRODUCTION 

At least ~ 20% of exoplanet s are associated with 
one or more stellar c ompan i ons dRaghavan et al.l 
Desidera fc Barbierl 120071: lEggenberger et al ' 



2006 



Harbienl iZUUYt lEggenberger et alJ 120071: 
Neuhauserl 120091 : iRaghavan fc Stepransl 



Mugrauer 

20101) . Stellar companions may significantly alter 



the planetary orbits around their partner on secular 
timcscalcs. Close-in giant planets t end to be found 
prefe r entially in binary stel lar systems (jZucker fc Mazehl 
120021: lUdrv fc Santosll2007t ). On the other hand recent 
studies suggest that the frequency of giant planets in 
close binaries (< 100 AU) is significantly lowe r then 
in th e overall population ([Eggenberger et ail 120081 
120111 ) . indicating that distant stellar perturbers may be 
important in the production of HJs. 

Recent measurements of the sky-projected angle be- 
tween the orbits of several HJs and the spins of 
their host stars have shown that mi salignment and 
even retrograde orbits ar e common (e.g. iGaudi fc Winnl 
I2007t iTriaud et all 120101) . If these planets migrated in 
from much larger distances t hrough their interaction 
with the protoplanetary d isk dLin fc Papaloizoul 119861 : 
iMasset fc Papaloizo u 2003), their orbits should have low 
eccentricities and inclinations (but see lLai et ail 120101 : 
iThies et al.l 120 111) . An alternative scenario that can ac- 
count for the retrograde orbits involves the secular in- 
teractio n between a planet and a binary stellar com- 
panion (IWu fc Murravl2"00llFabrvckv fc Tremaiii32007t 
IWu et al.H2007t iTakeda et alj|2008t iCorreia et al.ll2011l) . 

Many theoretical studies have investiga ted the role of 
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l20TTallbl : IVeras fc Tout! |20L1 iKratter fc Pereti |20H . 
Long-term stability requires that the system be suffi- 
ciently hierarchical. For an inner binary (stellar mass 
mi and planet mass 777,2) in a nearly-Kcplcrian orbit with 
semi-major axis (SMA) a\, the outer orbit for mass 7773 
around the center of mass of the inner binary must have 
SMA ci2 3> Q4 . For stability the eccentricity of the outer 
orbit, e2, must also be small enough for 7773 to avoid 
close approaches with the inner orbit. In such systems 
a sufficiently inclined 777,3 can produce large-amplitude 
oscillations of the eccentricities and inclinations of the 
777,1-777,2 system, while a\ and 02 remain nearly co nstant; 
this is the so-ca lled Kozai-Lidov mechanism (jKozail 

IT9^lLidovlll^6l . 

The standard treatment dKozail I1962D as sumes that 
7772 <IC mi, ™3 j and e2 ~ 0. Recently, INaoz et"ai1 
(|201 lal TbT) showed that these approximations are not ap- 
propriate for many systems. In the presence of an eccen- 
tric outer orbit, or if the test particle approximation for 
the inner planet is relaxed, the behavior of the system 
can be qualitatively different. The different behavior is 
associated with breaking of the axial symmetry present 
in the standard Kozai analysis; the associated lack of 
conservation of the axial component of the inner orbit's 
angular momentum allows the orbit to reach extremely 
high eccentricities and can even "flip" the orbit from pro- 
grade to retrograde (with respect to the total angular 
momentum) . 

Here we study the evolution of Jupiter-like planets in 
binaries through Monte Carlo simulations. We include an 
octupole-level approximation to the perturbing potential 
and the ti dal interactions between the star and t he planet 
(following lEggleton fc Kiseleva-Eggletonl 120011) , as well 
as GR precession. 

2. NUMERICAL SETUP 

We set mi = 1M , m 2 = IMj and 7773 = 1M 
for all our calculations. We denote the inclination an- 
gle of the inner (outer) orbit with respect to the total 
angular momentum by i\ (12), so that the mutual in- 
clination between the two orbits is i tot = ii + *2- We 
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name 


0-2 




ty 2 


HJs 


ip > 90° 


Total 




[AU] 


(initial) 


[yr] 


% 


% 


Runs 


SMA1000 


100 


0.01 


1.5 


13 


50 


2478 


SMA100 


1000 


0.01 


1.5 


1 


60 


3165 


SMA500 


500 


0.01 


1.5 


19 


11 


2130 


SMA500L 


500 


0.01 


0.015 


18 


16 


2014 


SMA500H 


500 


0.01 


150 


0.5 


25 


2514 


SMARan 


uniform in loj 


; 0.01 


1.5 


16 


11 


4690 


SMARanclR 


uniform in lof 


; Raylcigh 


1.5 


10 


39 


1674 



TABLE 1 

Parameters of the Monte Carlo simulations. In all runs we start the Jupiter in a circular orbit at 5AU. The arguments 

OF PERIASTRON, gi AND (?2, WERE CHOSEN FROM A UNIFORM DISTRIBUTION, AND THE MUTUAL INCLINATION BETWEEN THE INNER AND 
OUTER ORBITS IS DRAWN ISOTROPICALLY. On A FEW DAYS ORBIT, ty = 1.5 YR IS EQUIVALENT TO Q ~ 2 X 10 5 WHILE ty = 0.015 YR AND 
ty = 150YR (USED IN THE SMA500L AND SMA500H RUNS) ARE EQUIVALENT TO Q ~ 2 X 10 3 AND 2 X 10 7 , RESPECTIVELY. FOR THE STAR 

WE CHOOSE ty = 50 YR IN ALL RUNS, WHICH IS EQUIVALENT TO Q ~ 10 5 . 



call the angle between the spin of the inner star and 
the direction of the angular momentum of the inner or- 
bit ij). The projection of ip onto the plane of the sky, 
A, can be o bserved through the Rossiter-McLaughlin 
effect (e.g.. IGaudi fc Winn 120071) and other methods 



/e.g-INutzman et al.ll20lUISanchis-Oieda fc WinrJl2011t 
IShporer et al.ll2012t) . 

We solve the octupole-lev e l secula r equations numer- 
ically following iNaoz et"al"1 (l2011b[ ). We are able to 
follow the spin vectors of both the planet and the 
star. The effects of ti d es an d spins agree with those of 
iFabrvckv fe Tremainel (|2007f ) (private communication). 
In many of our simulations, the inner orbit reaches ex- 
tremely high eccentricities. During excursions to high 
eccentricity there is a competition between the increased 
efficiency of tides leading to the Kozai capture process 
(N aoz et aLll2011aD and the possibility of destroying the 
system by crossing the Roche limit. 

The upper limit for each system's integration time in 
all our simulations was 9 Gyr. However, if the planet 
becomes tidally captured by the inner star, the integra- 
tion becomes extremely expensive. We therefore stop 
the simulation, classifying the planet as a HJ, when- 
ever ai < 0.03 AU and ei < 0.01, or ai < 0.06 AU and 
ei < 0.4. These conditions imply either a circular HJ, 
or an eccentric one, respectively. If a planet crosses the 
Roche limit we also stop the run and we consider it lost 
(see below). 

We draw outer periods from the log-normal distribu- 
tion of IDuquennov fc Mayor! (11991D . but require that the 
corresponding a 2 lies between 51 and 1500 AU. At these 
separations the distribution is approximately constant 
in log(P 2 ), where P 2 is the outer orbit period. The 
distribution of the outer eccentricity is assumed ther- 
mal We also performed runs for which we fix a 2 = 100, 
500 or 1000 AU. These runs serve as a "zoom-in" since 
they give slices through the larger SMA parameter space. 
Two additional runs have different tidal friction param- 
ete rs, varied over two ord ers of magnitude. Motivated 
by iMoorhead et all (|2011[ ). we adopt a Rayleigh distri- 
bution for ei with a mean eccentricity of 0.175 for one 
of the runs. All runs are summarized in Table [T] Note 



4 Recently, Raghavan & Stcprans (2010) showed that the eccen- 
tricity distribution for sun-like stars is closer to uniform. Therefore 
we have conducted an additional run similar to SMA500 but with 
the stellar-binary eccentricity drown from a uniform distribution. 
We have found that the results were essentially unaffected. 



that all our initial conditio ns are stable according to the 
iMardling fc Aarsethl (|200l criterion. 

The differential equations that govern the in- 
ner binary's tidal evoluti o n we re presented in 
lEggleton fc Kiseleva-Eggletonl (|2001l ). These equa- 
tions take into account stellar distortion due to tides 
a nd rotation, with tidal dissipation based on the theory 
of lEggleton et al.l (|1998| ). The viscous time scale, ty, is 
related to the quality factor Q ([Goldreich fc Soterlll966D 
by 

n _ 4 k L Gm ty 

W ~ 3(l + 2fc L ) 2 B? n ' [ ' 

where n = 2tt/P is the mean motion of the or- 
bit and is the classical apsidal motion constant. 
We use the typical value = 0.014, valid for 

n = 3 polytropes, when representing stars and ki = 
0.25, valid for n = 1 polytropes, when representing 
gas giant planets (lEggleton fc Kiseleva-Eggletonl 120011 : 



IFabrvckv fc Tremaind 120071 ). Table HI specifies the dif- 
ferent viscous times we have used in our treatment for 
the Jupiter-like planet. 

The EKL mechanism can cause evolution of e\ to 
extremely high values, implying a high probability 
that a planet wi l l cros s the Roche limit. Following 
iMatsumura et all (|2010f ) we define 



Rl 



R2 
0.6 



m 2 



m 2 + mi 



-1/3 



(2) 



where i? 2 = lRj. If we find that ai(l — ei) < Rl we 
stop the run and assume that the planet is lost. 

3. RESULTS 

In Figure [T] we show the fraction of systems that have 
survived the EKL mechanism as a function of SMA (for 
both the SMARan and SMARanelR runs), and as a 
function of the outer orbit eccentricity in the "zoom-in" 
runs. In these figures we also show the fraction of sys- 
tems that formed HJs via the Kozai capture mechanism. 
For wide separations (a 2 > 100 AU) the EKL mechanism 
and Kozai capture are far more efficient than suggeste d 
previously (|Wu et al.ll2007l : IFabrvckv fc Tremainell2007D . 
In Figu re [H left panel, we also compare our results to 
those of INaoz et al.l (|2011aD for a planetary perturber. 

In Figure [2] we show the initial conditions that are asso- 
ciated with HJ formation. Most of the systems that form 
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Fig. 1. — Fraction of planets as a function of outer SMA (left) and the outer orbit eccentricity (right). In the left panel we consider the 
SMARan (solid lines) and SMARanclR (dashed lines) runs. We show the fraction of systems that survive the EKL process (blue lines), the 
fraction of systems that formed HJs via Kozai capture (red lines) and the fraction of HJs that ended up in retrograde motion with respect 
to the stellar spin axis (ip > 90°, brown lines). In the bottom right panel we consider the SMA500 run (solid thick line), the SMA1000 
run (solid thin line), the SMA100 run (thin dashed lines), and the run with extremely efficient tides, SMA500L (dot-dashed lines). The 
SMA500H run had extremely low efficiency of HJ formation and does not appear here. In t he top right panel w e compare the SMA500 run 
(solid thick line) to a system in which the perturber is another giant planet (dashed lines) (Naoz ct al. 2011a). 




Fig. 2. — Formation of HJs as a function of initial conditions. We 
consider the initial mutual inclination, itot,0i an d the initial outer 
orbit eccentricity, e2,o- The complete suite of runs are marked by 
the light blue dots. In dark blue we mark the systems that form 
HJs, and in red we show the systems that did not survive the EKL 
process. The solid line represents the analytical prediction for the 
first flip (see text). The different runs we considered are labeled in 
the figure. 



HJs lie within the bounds of an analytic criterion for de- 
termining the regions of parameter space associated with 
flipping the inner orbit in the test particl e approximation 
(|Lithwick fc NaozteOlWKatz et al.l201lD . Before the flip 
occurs, the inner orbit's eccentricity becomes extremely 
high, often resulting in Kozai capture. Overall, the EKL 
mechanism and Kozai capture produce HJs about 15% 
of the time for the systems we considered. 

In Figure [3] we show the distribution of mutual inclina- 
tion, i to t) and spin-orbit angle, tp, resulting from our sim- 
ulations. Consideration of the octupole-level secular evo- 
lution when the inner body is a test particle (|Katz et al.l 
1201 ID shows that the average part of the octupole Hamil- 
tonian peaks for minimum inner eccentricity of ~ 0.335 
and inclination i to t — 61.7°. This may be related to 
the peak observed in Figure [3] around 60° . An identical 
peak occurs for the test particle octupolar Hamiltonian 
at the corresponding retrograde inclination of 120°. In 
this work we find that the retrograde peak is wider, sug- 



gesting that r etrograde orbit s are more stable than pro- 
grade orbits (jlnnanenl 119791 119801: iMorais fc Giuppond 
l20ll . 

4. SPIN-ORBIT ANGLES AND COMPARISON WITH 
OBSERVATIONS 

Any proposed mechanism for producing misaligned 
HJs should be consistent with both the observed shape of 
the spin-orbit angle distribution and the overall rate of 
HJ formation. Here we examine these questions for HJs 
produced via the EKL mechanism and Kozai capture. 

The efficiency of HJ formation from the EKL mecha- 
nism for close (d2 < 500 AU) binary systems is lower for 
more distant perturbers, in accordan ce with observations 
(e.g., lEggenberger et al.l 120081 120111) . The formation ef- 
ficiency for wide binaries is almost independent of both 
SMA and eccentricity for t he outer orbit (see Figure [lj . 
Following IWu et al.l (pOOl we estimate the fraction of 
stars with HJs as 



(3) 



where fb is the fraction of stars in binary systems, f p 
is the fraction of stars with Jupiter-mass planets formed 
at a few AU, and /ekl — 0.15 is the efficiency of our 
mechanism from simulations. Taki ng ft, ~ 0.3, estimated 
from t he tail of the distribution in lDuauennov fc Mayor! 
(119911) . and f„ ~ .07, consistent with the studies from 
iWright et all 120121 ) and lMarcv etall ([2001 . we find / ~ 
0.0032. The fraction of stars hosting HJs is about 1%, 
so we estimate that our mechanism can provide ~ 30% 
of the total HJ formation. 

Figure 2] shows the observed projected spin-orbit angle 
distribution, the distribution produced by our mecha- 
nism, and the distribution produced by the dynamical 
planet -planet scattering considered in iNagasawa & Ida 
(|201ir FL It is clear that no one mechanism is a 
good fit to the observed distribution: more systems 
are observed to be aligned than are produced by any 



5 See also the recent c alculations bv IBolev et al. (20l2) and 
IBeauge fc Nesvornvl pHlT ). 
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Fig. 3. — Effect of initial conditions on the final state for HJs formed in the SMARan run. In panel (a) we show the distribution of the 
initial mutual inclination (itot) for the systems that formed HJs. We also show in light colors the initial mutual inclination for those HJs 
that ended up in retrograde spin-orbit configuration (i.e., ip > 90°). In panel (b) we show a scatter-plot of the initial, ii,o, and final, 11,1? , 
inner inclinations for the SMARan run. Panel (c) shows the distribution of the final mutual inclination for the SMARan run (thick solid 
blue line); we also show the distribution of the final mutual inclination for the SMA500 and SMA1000 runs (thin solid red and dashed black 
lines, respectively) as a function of the initial «tot,o- Panel (d) shows the distribution of the final spin-orbit angle, i/>, for the SMARan run 
(thick solid blue line). We also show the distribution of the final ip for the SMA500, SMA1000 and SMA500L runs (thin solid red, dashed 
black and dot-dashed green lines, respectively). 



of the proposed formation mechanisms. Aligned sys- 
tems ar e a more natura l con s equence of disk mi- 
grati on dLin fc Papaloizoul 119861: iM assct & Papaloizoul 
120031) . Following iMorton fc Johnson! (|2011[ ). we express 
the complete spin-orbit angle distribution as a sum of 
contributions from an aligned component, an EKL com- 
ponent and a dynamical plane t-planet scattering com- 
ponent (jNagasawa fc Idall20Tl 

P(M {fi}) = (1 - /EKL - /scattering) S(X) (4) 
+ /ekl PEKl(^) /scattering ^scattering (A) , 

with 5 the Dirac delta function and pekl(A) and 
Pscattering(A) the distributions shown in Figure |4] The 
fi represent the relative contribution of each component 
to the total distribution. We impose uniform priors on 
the fi with fi > and J^. /, < 1, and compute the poste- 
rior probability of the set of fi, {fi}, given the observed 
spin-orbit alignments (and associated errors) of 45 sys- 
temt0 from lWright~et~aTl (|20TlT) . d, using Bayes rule: 



P ({/*} |d) «p(d|{/i})p ({/<}). 



(5) 



Treating each spin-orbit angle measur e ment as a Gaus- 
sian with errors given bv lWright et"aLl (f20TTh . the likeli- 
hood for the data set, d, is given by 



45 „ 

P(d\{fi}) = l[ I dXN (A; A. i; of) p (A {fi}) , (6) 

i— 1 

where the Xi are the observed angles, of the observa- 
tional errors, and N(x; [i, a) is the Gaussian PDF with 



6 The spin-or bit observations reported at 

http : //exopla n ets . o rg (Wright ct al. 2011] ), include da t a from 
Simpson ct al. (2011); Narita ct al. (2010); Triaud ct al. (2010); 
Queloz ct al. (2010); Gucntlicr ct al. (2012); Triaud cTaLl ip009ft : 
Winn et alT (120091); I Johnson et all (120091) ; I Jenkins et all COW); 
MoutouetaLI (120111); [Wolf et all (120071); IGillon et alT EuOa) ; 
Winn et al.l (12011 . 2005): Brown ct al. (2012): Whm etaD (120071. 



20101) ; IJohnson et all (2008); Hcllicr ct al. (2011); Bouc hv et al 
( 20081) ; INarita et all (12 009): Hc brard etHTl ( 12011 ): Simps on et al 
( 2rjTol) : Anderson ct al. ( 20111): INarita et all p007l '): IHehrard et al 
(2010): Pont ct al. (2010). 



mean p and standard deviation o evaluated at x. The 
integral is evaluated over the range A G [O.tt], properly 
accounting for wrapping at the endpoints. The result- 
ing probability distribution^ for the fractional contribu- 
tion of each component (aligned, EKL, and planet-planet 
scattering) appear in Figure |U We also show in Figured] 
the probability distribution for the EKL contribution to 
the misaligned systems, 



^misaligned 

JEKL 



EKL 



EKL 



scattering 



(7) 



As expected from the distributions in Figure |U the 
data support a significant aligned component, account- 
ing for about 50% of the observed systems. The posterior 
is nearly degenerate along the line /scattering + /ekl — 0.5 
corresponding to a total contribution from scattering and 
EKL of 50%. However, the data do prefer a larger con- 
tribution from EKL than scattering, with EKL account- 
ing for most likely ~ 30% of the observed systems and 
planet-planet scattering ~ 10% to 20% of the systems. 
The data prefer that EKL produces between 60% and 
80% of the misaligned systems, but fractions as low as 
0% and as high as 100% cannot be ruled out. Note that 
the allowed contribution from EKL to the shape of the 
spin-orbit angle distribution is consistent with the es- 
timate from the rate of HJ production that EKL can 
account for about 3 0% of all HJ systems. 

Previous studies (fFabrvck v fc Trema inc 2007]) of the 
secular effects of stellar perturbers have considered only 
the quadrupole terms in the secular potential. The ex- 
pansion up to octupole order used in our work pro- 
duces qualitatively different behavior. The additional 
terms in the potential can drive the inner orbit to 
much more extreme eccentricities and inclinations (in- 
cluding retrograde inclinations, which are impossible in 
the quadrupole limit), and, based on Figure [1] leads 

7 This analysis was performe d using customized code implement- 
ing the algorithm described in Porcman-Mackey ct al. (2012) and 
references therein. 
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Fig. 4. — (Upper Left) Distributions of proj ected (i.e., observed) spin— orbit angle, A: in black, the observed distribution of spin— orbit 
angles in the 45 systems of j Wright ct al. ( 2 0111 ') (the distribution has been smoothed using a Gaussian kernel that appropriately takes into 
account the observational errors in each system) ; in green, the distribution of spin-orbit angles produced by the EK L mechanism discussed i n 
this work; and in red, the distribution of spin-orbit angles produced by the planet— planet scattering simulations inlNagasawa fc Ida! 1120111 ), 
(Upper Right) Posterior distributions of the fractional contribution of the various components to the overall spin-orbit angle distribution 
implied by the data (see eq. and 10): in green, distribution of the fraction of the spin-orbit distribution due to the EKL mechanism 
(/ekl) implied by the data; and in red the distribution of the fraction due to planet-planet scattering (/scattering : INagas awa & Idal lMtTr) ; 
and in blue the fraction attributed to disk migration (1 — /ekl — /scattering)- (Lower Left) Posterior probability distribution in the 
/scattering _ /EKL plane. (Lower Right) Probability distribution for the fraction of misaligned systems produced by EKL (sec cq. 171). 



to much more efficient HJ formation than found in 
iFabrvckv fc Tremaind (|2007t) . 

5. SUMMARY AND DISCUSSION 

We studied the formation and evolution of HJs in wide 
stellar binaries using an approximation accurate to oc- 
tupolc order in the SMA ratio. Recent studies have 
shown that octuple-level perturbations can play a very 
important role in t he dynamics of thre e-body systems 
(jNaoz et all 1201 lbl ). INaoz et all (|2011aD showed that in 
the presence of a second planetary-mass perturber, sec- 
ular perturbations can easily produce retrograde orbits. 

From an observational point of view, statistical analy- 
ses suggest that at least ~ 20% of the known extrasolar 
planetary sys tems are associated with one or more stellar 
companions (iRaghavan et al.l 120061 : iDesidera fe Barbieril 
120071 : lEggenberger et al.l l2007t ). and thus stellar per- 
turbers are likely to be at least as important as planetary 
perturbers. 

Our results differ from those of previo us studies 
(|Fabrvckv fc Tremaind [20071 IWu et aljf2007l ) conducted 
at quadrupole order. We have found that the EKL mech- 



anism produces HJs through Kozai capture about 15% of 
the time. Given realistic assumptions about the rate of 
binary stellar systems and the fraction of systems that 
initially host planets, we find that our mechanism can 
produce HJs in about 0.3% of stars. Since about 1% of 
stars host a HJ, EKL may account for about 30% of all 
HJs. 

By comparing the shape of the sky-projected spin- 
orbit angle distribution p roduced by EKL, dyn amical 
planet-planet scattering (|Nagasawa fc Idal l20ll . and 
disk migration, to the observed spin-orbit angle distribu- 
tion, we find that EKL likely contributes to about 30% 
of the observed distribution, consistent with the overall 
rate of HJ production. 
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